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We calculate the single-particle spectral function for the one-band Bose-Hubbard model within 
the random phase approximation (RPA). In the strongly correlated superfluid, in addition to the 
gapless phonon excitations, we find extra gapped modes which become particularly relevant near 
the superfiuid-Mott quantum phase transition (QPT). The strength in one of the gapped modes, 
a precursor of the Mott phase, grows as the QPT is approached and evolves into a hole (particle) 
excitation in the Mott insulator depending on whether the chemical potential is above (below) 
the tip of the lobe. The sound velocity c of the Goldstone modes remains finite when the transition 
is approached at a constant density, otherwise, it vanishes at the transition. It agrees well with 
Bogoliubov theory except close to the transition. We also calculate the spatial correlations for 
bosons in an inhomogeneous trapping potential creating alternating shells of Mott insulator and 
superfiuid. Finally, we discuss the capability of the RPA approximation to correctly account for 
quantum fiuctuations in the vicinity of the QPT. 

PACS numbers: 05.30.Jp, 03.75.Lm, 71.45.Gm, 37.10. Jk 



I. INTRODUCTION 

Optical lattices have made it possible to explore the 
properties of ultracold dilute atoms in a new regime of 
strong correlationsii^i^. By tuning the strength of the 
laser field, the effective interactions between atoms can 
be tuned to become stronger than their kinetic energy. 
For Bose systems, such a competition between kinetic t 
and interaction energy U drives a quantum phase transi- 
tion"*'^ from a kinetic energy dominated superfluid (SF) 
phase to an interaction dominated Mott insulating (MI) 
phase. The Bose Hubbard model (BHM) captures the 
essential physics of this problem^, provided the interac- 
tions between the bosons are smaller than interband en- 
ergies and the problem can be treated in the single band 
approximation. While the BHM was proposed much be- 
fore optical lattice experiments became available, a direct 
experimental realization was missing. In condensed mat- 
ter systems, Josephson junction arrays^, ^He in vycor 
and aerogels^, vortices in superconductors* and quan- 
tum magnets^ can be modeled by the BHM, but the ac- 
tual systems have additional complications of disorder 
or longer range interactions which make the comparisons 
between theory and experiment difficult. 

One of the main advantages of the cold atom systems 
is that they are clean and much more tunable: the den- 
sity of bosons, their effective interaction, the tunneling 
amplitude between the wells, the number of lattice sites, 
the shape of the trapping potential and aspect ratios can 
all be varied rather easily, making it possible to study 
the effects of inhomogeneity and dimensionality. In ad- 
dition, it is possible to add random potentials to study 
the effects of disorder. This sets the optical lattice sys- 
tems apart as a useful testing ground for theoretical ideas 
in the area of strongly interacting bosons and fermions. 
This model has also provided tremendous impetus for the 
development of new measurement techniques to address 



questions about the nature of the excitations especially 
near the transition 22AL12J^J^^1^J^JJ^. The recent ex- 
periments on the dynamic o^^i^° have given a window into 
the different time scales operating within the different 
phases and around the quantum phase transitions. 

The paper is organised as follows: In Sect. HH we 
present the Bose-Hubbard model and state our main re- 
sults. In Sect, mil we discuss the nature of the spec- 
tral function calculated within the RPA formalism as 
it evolves from the SF to the MI phase by decreasing 
t/U. The sound velocity in the SF phase and its com- 
parison with Bogoliubov theory is contained in Sect. IIVI 
The momentum distribution and the spatial correlations 
are calculated in Sect. |Vl The RPA formalism is gen- 
eralized to a spatially inhomogeneous trapping potential 
in Sect. IVII We conclude in Sect. IVIII with some re- 
marks about the comparison between RPA and mean 
field theory. There are three appendices that give the 
details of the calculations of the Green function within 
RPA (App. E| . the Bogoliubov calculation for the Bose- 
Hubbard model (App. |B]), and the momentum distribu- 
tion function within RPA in the Mott regime (App. [C)) . 

II. MODEL AND MAIN RESULTS 

The Bose-Hubbard Hamiltonian is defined as 

H = -^ Yl ("^"J +a,a])+^^ni(nj-l)-/i^n, , 

(ij) i i 

(1) 

where and a] are bosonic annihilation and creation 
operators respectively and Ui — a|ai is the density op- 
erator. The parameter U describes the on-site repulsive 
interaction between bosons, t is the tunneling parameter 
between nearest neighbors as indicated by the symbol 
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(ij) , fi is the chemical potential that fixes the number of 
particles and z — 2D is the coordination number in D 
dimensions. This Hamiltonian shows a quantum phase 
transition (QPT) from the SF to the MI phase as a func- 
tion of t/U. The theoretical approaches used to investi- 
gate this Hamiltonian include mean field theory^l^ per- 
turbation theory^^, variational methods^"^ and quantum 
Monte Carlo simulations 24.25^ 

In this paper, we use a Green's function formalism to 
investigate the excitations and correlations in the BHM. 
We start from the mean-field (MF) ground state which is 
essentially a product of single site states and go beyond 
it by including inter- well coupling within a random phase 
approximation (RPA)— i^. 

Our main results are as follows: 

(1) In the weakly interacting SF {t/U « 100) the gap- 
less phonons of the SF, or the Goldstone modes arising 
due to the broken gauge symmetry, exhaust the sum rule 
on the total spectral weight, as expected. For t/U ~ 10 
there are already small deviations from Bogoliubov the- 
ory and new gapped modes appear in the SF phase. 
These gapped modes pick up strength ds, t/U « 1. The 
sum rule is now satisfied only upon including both phonon 
and gapped modes. 

(2) At the transition, we observe the progression of one 
of the phonon modes in the SF into a gapped mode in 
the MI (the one which is gapless at the QPT). The sec- 
ond gapped mode in the MI instead arises from one of 
the gapped modes in the SF. Such gapped modes in the 
SF have been reported previously using several theoret- 
ical methods I4^i5,i6^i7 _ ^j-g^g ^jj^t these additional 
gapped modes are a distinctive signature of a strongly 
correlated SF in proximity to a MI in an optical lattice. 
They indicate the redistribution of spectral weight from 
the coherent phonon modes into incoherent excitations, 
and are a precursor of the MI beyond the QPT. 

(3) We calculate the sound velocity in the RPA for- 
malism and show that it agrees with c — 1/^ Km* cal- 
culated independently from the mean-field effective mass 
and compressibility. In a wide range of parameters, ex- 
cept very close to the SF-MI phase transition, the above 
sound velocity compares well with the predictions of Bo- 
goliubov theory. 

(4) We exploit a special feature of superfluids that al- 
lows us to extract the condensate fraction no from the 
strength of the phonon modes in the spectral function. 

(5) We calculate the spatial correlations in the case 
where an inhomogeneous confining potential is superim- 
posed on the optical lattice. The response to a perturba- 
tion is strongly influenced by the presence of alternating 
shells of Mott insulator and superfluid regions. 



III. FORMALISM: RPA APPROXIMATION, 
SPECTRAL FUNCTION, AND EXCITATIONS 

We start with the mean-field approximation in real 
space^i obtained by giving the annihilation and creation 



operators an expectation value defined by (a) = {a^) = ip. 
The order parameter tp identifies the nature of the sys- 
tem: it is non-zero in the SF phase and vanishes in the 
insulating phase. Substituting a — ip + = ip + a\ 

in Eq. ([1]), where a and are the fluctuations of the 
Bose field around the mean field value, the Hamiltonian 
H can be rewritten as a sum of on-site Hamiltonians 



Hf''^ = -^n,{n, - 1) - fin, - tip{al + a,) + t^^ , (2) 

which include the tunneling at the mean field level 
through the order parameter tp. In the MF approxima- 
tion, wc neglect the non-local inter-well hopping term 

~{t/2z) (^o-lcij + flifl]^ , which wc will later treat in 

RPA. 

The H^ ^^ can be diagonalized numerically, leading to 
a set of on-site eigenstates such that Hi\ia) = ecjia). 
In the Mott limit the eigenstates \ia) are number states, 
while in the SF regime they are coherent superpositions 
of several number states, allowing the order parameter 
to be different from zero. The MF ground state solution 
is given by the product state |$) = H^ji, 0), equivalent to 
the one obtained in the Gutzwiller Ansatz, where |z, 0) is 
the ground state of . 
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FIG. 1: (Color online) Mean field phase diagram in the ^/U 
vs. t/U plane. The blue line shows the Mott insulating lobe 
at density n = 1. On this diagram, we indicate two points 
where the QPT happens, that we are going to discuss in this 
paper: a generic one at n/U = 0.5 and t/U « 0.167 (black) 
and the tip of the lobe at fi/U = ^2 - 1 and t/U ^ 0.1716 
(red) where the QPT takes place at constant density. 



Within the mean field approximation, the state of the 
system is described by a product state over the differ- 
ent wells, neglecting all inter-well correlations. However, 
even in the MI, correlations between neighboring wells do 
not vanish; in fact they get large as t/U is increased from 
the Mott side, ultimately diverging at the transition. Ex- 
perimental evidence of these correlations is found in the 
interference picture of an atomic cloud released from a 3D 
optical lattice, the visibility of which does not suddenly 
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vanish at the phase transitioii22ii22ii2£i^i2^. These impor- 
tant features are captured by the RPA performed on the 
non-local tunneling terms of the BHM. At the end of this 
paper, we discuss the limitations of the RPA method and 
how it compares with the mean-field approximation. 

To go beyond the mean-field approximation, we treat 
the inter-well hopping term within RPA, as described 
in App. I^^'^^. This method allows us to compute 
the Green's function G(q, w) = {{a^^; Qq)) u , defined 
in Eqs. (|A4IA6p . and from that, the spectral function 
^(q,c^) = -(l/7r)ImG(q,c^). 

Due to the commutation relations of the bosonic de- 
struction and creation operators, the spectral function 
always satisfies the sum rule 



^(q, Uj)dLj — 1. 



(3) 



From the spectral function, one can extract the excitation 
spectrum, the strength of the excitation modes, and the 
related density of states 



DOS(w) 



^(q, uj)dc[. 



(4) 



Moreover, the spectral function is an essential ingredient 
to compute the momentum distribution 



A. Deep Mott regime 

In the deep Mott regime (zero tunneling) , for U {n 
1) < /i < Un, one finds 
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t^M7(q) 

riAiiiq) 



uj — {Un — fi) Lu — (U {n — 1) — iJ.) 

(7) 

Un~'fi>0 
[/(n-l)-^<0 

n, V q, (9) 



where n is the atomic occupation at each lattice site. The 
spectral function consists of two 5-functions, one at pos- 
itive energy Un — (relative to the chemical potential) 
required to add a particle and one at negative energy 
U{n — 1) — /i, required to remove a particle or add a 
hole to the MI, as seen in Eq. ([5]). The spectral function 
A{c[,uj) obtained using expression ([7]) trivially satisfies 
the sum rule in Eq. ([3]). The momentum distribution, as 
defined in Eqs. (|5l9p . is completely flat, corresponding to 
vanishing sitc-to-site correlations, and normalized to the 
total number of atoms in the lattice (n times the number 
of sites). Correspondingly, the single particle density ma- 
trix ^ , given by the Fourier transform of the momentum 
distribution, shows strictly on-site correlations. 



n(q) = (aljOq) = / y^(q,w)da 

J — OO 



(5) 



and the single particle density matrix given by the 
Fourier transform of the momentum distribution in real 
space 



P(r,r') = 



q 



'n(q), 



(6) 



where N indicates the number of lattice wells . The long 
distance behavior of p(r, r') as a function of the relative 
distance approaches the condensate density no which is 
non zero in a SF and vanishes in the MI. In the following 
we will calculate and discuss all these quantities. 

A fundamental implication of broken symmetry for 
bosonic systems is that the Goldstone modes are directly 
reflected in the single particle spectrum. In other words, 
phonons which are related to modes of density-density 
fluctuations (or two-particle Green function) also show 
up as the poles in the single particle Green function^. 
We study the behavior of the poles of the Green's func- 
tion, their strength, momentum and frequency depen- 
dence, to extract information about the excitations of 
the system. In the two extreme limits of deep MI and 
weakly interacting (Bogoliubov) SF, the Green's function 
can be calculated analytically and from it, the excitations 
frequencies and the momentum distribution r7,(q). 



B. Weakly interacting regime 

In the weakly interacting SF regime^, we have 



Gi3G(q,w) 
WBG(q) 

nBG(q) 



1 



(10) 

(11) 
(12) 



where Uq and Vq are the Bogoliubov amplitudes, LUq is 
the Bogoliubov frequency at momentum q (see App. [B] 
for details), and no the condensate density. In the 
weakly interacting SF regime, the sum rule in Eq. ([3]) 
is constrained by the Bogoliubov normalization condi- 



tion ImoP - 



1. The excitation energies are given 



by symmetric poles at positive and negative frequencies 
corresponding to the energies of the Bogoliubov spec- 
trum, highlighting in particular phononic excitations at 
low momentum. 

The momentum distribution, given by integrating the 
spectral function over negative energies, has a singular 
contribution from the condensate at zero momentum. 
The integral over all momenta different from zero gives 
the number of non condensed atoms, contributing the de- 
pletion from the condensate. In the regime where Bogoli- 
ubov theory is valid, the depletion is negligible compared 
to the atoms in the condensate. 
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C. Progression from SF to MI 

The RPA formalism allows us to calculate the spec- 
tral function with special emphasis on the strong cor- 
relation region near the QPT. In the deep SF, we find 
phonon collective modes reflected in the single particle 
spectrum. As t/U is decreased, the spectral weight is 
redistributed over a multi-mode structure composed by 
coherent phonon excitations and gapped single-particle 
excitations. When entering the MI phase at the QPT, 
the spectral weight reorganises and is shared by only two 
gapped modes, describing single particle excitations, one 
at positive and one at negative energy. In the following, 
we will deiscuss this behaviour more in detail. 



We use the position of the poles of the Green's function 
to determine the following results about the excitations 
of the system in the different regimes: 

(i) For a large number of particles per site (« 100) 
and weak interactions [t/U ~ 100), we exactly recover 
the Bogoliubov results. We point out that being able to 
describe the weakly interacting regime starting from the 
BHM is not a trivial result, because of the large number 
of basis states required (almost 150 states per site). 

(ii) By increasing the interactions and decreasing the 
number of particles per site, we observe small deviations 
between the spectrum obtained by RPA and the Bogoli- 
ubov theory: additional modes appear at higher frequen- 
cies, as shown in Fig. HJa) ior t/U = 10, in contrast 
with Bogoliubov theory which predicts a single excitation 
mode. While there are also differences in the dispersion 
of the sound modes at large momenta, we find in general 
that the Bogoliubov prediction turns out to be quite ac- 
curate in describing the low-q part of the spectrum and 
the sound velocity, even in the case of strong interactions 
(see Sect. HV)) . 

(iii) As t/U becomes of order unity and the effect of 
strong correlations grows, additional gapped modes in 
the SF phase are clearly visible and grow in strength as 
seen in Fig. [^^b). The phonon modes are not sufficient 
to exhaust the sum rule in Eq. In a strongly inter- 
acting SF (e.g for t/U < 0.25, as shown in Fig. [DJc)), 
many modes (in the cases we have considered, up to four 
at positive and four at negative energy) have to be in- 
cluded in order to exhaust the sum rule in Eq. ([3]). In 
the particular case of /i/C/ = 0.5 and t/U = 0.25, the full 
dispersion of the modes is shown in Fig.[3l^b). 

(iv) In the MI , only two excitation modes exist, as 
shown in Figs.^d) and[3]^a). The mode at positive and 
the one at negative energy correspond respectively to the 
energy needed to create a particle or a hole in the sys- 
tem. For any given t/U, the difference of the excitation 
energies at g = exactly coincides to the width of the 
mean-field Mott lobe at the same t/U. 
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FIG. 2: (Color online) Spectral function A{q, ui) as a function 
of ij for various q. Results obtained by RPA (black) and 
Bogoliubov theory (green), (a) Weakly interacting SF: t/U — 
10 with ~ 10 bosons per site; the RPA calculation agrees 
extremely well with the Bogoliubov theory; notice indications 
of additional modes at higher ui. (b) t/U — 1 with « 1.8 
bosons per site, (c) Strongly interacting SF: t/U — 0.25 with 
~ 1.1 bosons per site; stronger deviations from Bogoliubov 
theory are present especially at larger q. Additional modes 
are clearly visible in the spectrum, (d) Mott insulating phase 
for t/U = 0.16 and 1 boson per site. In all those figures 
fi/U = 0.5. 



D. Strengths of the spectral function 



The progression of the modes from the strongly corre- 
lated SF into the MI is better understood by calculating 
the strengths of the excitations Si, defined as follows: 
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FIG. 3: (Color online) Dispersion (a,b) and strength (c,d) 
of excitation modes at ^/U = 0.5: (a) and (c) are in the 
Mott regime t/U — 0.1; (b) and (d) are in the SF regime 
t/U = 0.25. For clarity, in (b) the 4th pair of resonances at 
u ~ ±18.4 is not shown and in (d) only the strength of the 4 
modes at lower energy is shown. Note that modes at positive 
(negative) energy have positive (negative) strength. 
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FIG. 4: (Color online) Energy (a) and strength (b) of the 
modes at low q = 7r/50 as a function of t/U for fj,/U = 0.5. 
Notice the presence of both phonons and gapped modes in the 
strongly interacting SF and their evolution into two gapped 
modes into the MI. One of the gapped modes in the MI evolves 
from the phonon mode in the SF and the other one from a 
gapped mode in the SF. The thin vertical line at t/U ~ 0.167 
indicates the QPT. 



A{q,Lo) = J2S^SiLU-LU^). (13) 

i 

Numerically, a small but finite imaginary part of the en- 
ergy regularizes the spectral function and provides an ac- 
curate fitting procedure to determine the position of the 
poles and their strength. We checked that the sum rule 
in Eq. which using Eq. implies J2i Si = 1, was 
found to be satisfied to better than few parts in 10~^ 
for all t/U. We are therefore confident that we have 
identified all the excitations which contribute in a non- 
negligible way to the spectrum. 



In Figs. 14151 we plot the position of the resonances and 
their strengths for a fixed value oi q = 7r/50, varying 
the parameter t/U across the phase transition for fixed 
chemical potential. As explained above, the many-mode 
spectrum in the SF phase evolves into the two-mode ex- 
citation spectrum in the MI. 

The transition from the MI to the SF phase occurs 
when one of the two Mott branches becomes gapless (see 
Fig. m^a)). This is the particle (hole) gapped mode in 
the MI depending on whether the chemical potential fi is 



above (below) the tip of the lobe. This Mott branch 
evolves into the phononic mode at positive (negative) 
energy, while the other phononic mode arises at zero 
strength without having a precursor in the Mott phase. 
The second Mott branch evolves into a gapped superfluid 
mode, and symmetrical in energy a second superfluid 
gapped mode arises with zero strength without having 
a precursor in the Mott phase (see Fig. [^b)). 

The behavior at the tip of the lobe is quite interesting. 
In that case both Mott gapped modes become simultane- 
ously gapless at the QPT (see Fig.[SJa)). From them the 
lowest four modes in the SF arise, two of them becom- 
ing phononic modes and two of them becoming gapped 
modes when moving away from the transition. However, 
within our approach, we flnd a similar behavior to the 
one described above, namely that one of the Mott modes 
evolves into a SF phononic mode, while the second one 
evolves into a SF gapped mode (see Figs.l^b)). 

A similar result was found by Huber and collabora- 
tors -ii using an effective 3-state approximation and map- 
ping onto a spin-1 Hamiltonian. They further suggested 
that measurements of the dynamical structure factor us- 
ing Bragg spectroscopy and lattice modulations could be 
an effective way to investigate the different modes. 

The problem was also investigated by Sengupta and 
Dupuis in the strong coupling regime by deriving an 
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FIG. 5: (Color online) Energy (a) and strength (b) of the 
modes at low q — 7r/50 as a function of t/U at fj^/U = %/2 — 1 
corresponding to the tip of the lobe. The frequency of the 
lowest four modes (two phononic and two gapped) in the SF 
vanish at the QPT. The thin vertical line at t/U ^ 0.1716 
indicates the QPT. In the inset, the zoom around the QPT 
of panel (b) is shown. 



effective action using Hubbard Stratonovich transforma- 
tions. By expanding the action to quadratic order in the 
fluctuations they found gapped excitations in the MI and 
gapless Goldstone modes in the SF. They found two ad- 
ditional gapped modes in the SF, which present a similar 
behavior to the one discussed in this paper. 
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E. Density of states 

A further quantity that one can use to characterise the 
excitations of the system is the density of states defined 
in Eq. (|4])^ii^. We calculate it across the QPT for a 
ID system^, as shown in Fig. [6l In Fig. [6lja) one can 
recognize the multi-mode structure of a strongly corre- 
lated SF through a clearly enhanced DOS in the energy 
range of the corresponding excitation branch. In par- 
ticular, we can see here two phononic branches and two 
gapped ones. When approaching the QPT {t/U — 0.17 
and Fig. EJJb)), we encounter a similar structure, where 
the width of the gapped branch at negative energy is 
increased, while the gapped branch at positive energy is 
not visible on the scale of this picture (although existing), 
since its strength goes to zero at the QPT. As expected, 
in the Mott regime, the DOS is different from zero only in 
the energy range of the two gapped excitation branches, 
one at negative and one at positive energy (Figs. [ni^c,d)). 



FIG. 6: Density of states for several values of t/U. (a) t/U — 
0.25, SF regime; (b) t/U = 0.17, SF regime, very close to the 
QPT; (c) t/U = 0.16 MI regime, very close to the QPT; (d) 
t/U = 0.1 MI regime. All those results are for fi/U — 0.5. 

As observed in Fig. EJc), the branch at positive energy 
extends almost to = 0, indicating the disappearance of 
the gap at the QPT. 



IV. SOUND VELOCITY 

The presence of phononic modes in the excitation spec- 
trum is an important signature of superfluidity. These 
modes disappear in the Mott phase, where sound cannot 
propagate because of a gap in the spectrum. In this sec- 
tion, we discuss the evolution of the sound velocity in the 
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strongly correlated SF phase as the SF-MI transition is 
approached. 

The sound velocity is related to the compressibility k 
and the effective mass j77, *36, 37^38 through the relation 



(14) 



where = p{dp/dp) and the SF fraction Ps/ p = 

m/m* . We calculate the sound velocity c from the slope 
of the gapless mode in the RPA spectrum 




lim Lij{q) 

q| ^0 



c|q|: 



(15) 



and compare it with the one obtained with the compress- 
ibility relation in Eq. ^4]) . We find perfect agreement 
between the values for the sound velocity extracted by 
the two different methods. It is important to note that 
the method using the compressibility relation in Eq. (fH)) 
only requires the knowledge of the mean-field solution, 
which provides the equation of state /i(p) and the SF 
density ps = 

We find that when the SF-MI transition, tuned hy t/U 
and fi/U, is approached at a generic point away from the 
tip of the lobe, the sound velocity vanishes, as shown 
in Fig. [TJa). This is due to the fact that at the transi- 
tion the compressibility remains finite, but the SF den- 
sity vanishes. Instead, the tip of the lobe where the phase 
transition happens at constant density and dp/ dpi, = is 
a special point: there, a perfect compensation between 
the divergent inverse compressibility and the vanishing 
SF density takes place, which results in a finite sound 
velocity as seen in Fig.[71[b). 

We complete our analysis by comparing the sound ve- 
locity calculated above with the results of Bogoliubov 
theory, which for a tunneling parameter t, on-site inter- 
action U and coordination number z, predicts the value 



2t 

CBG 



(16) 



as explained in detail in App. |B] The Bogoliubov predic- 
tions are remarkably good in a wide range of parameters 
and fail only in the close proximity of the phase transi- 
tion (see thin lines in Fig. [7]), since Bogohubov theory 
does not account correctly for the vanishing of the order 
parameter at that point. 



V. MOMENTUM DISTRIBUTION AND 
SPATIAL CORRELATIONS 

From Eq. ([5]), the momentum distribution n(q) is ob- 
tained by integrating the spectral function over negative 




FIG. 7: (Color online) Sound velocity c (blue full line) ex- 
tracted from the slope of the dispersion of the phonon modes 
and from Eq. 1141 order parameter (red dashed line) . Also 
shown for comparison are the sound velocity and the order 
parameter obtained from Bogoliubov theory (thin lines). In 
the inset dp/djj. is shown, (a) jj./U = 0.5; (b) p/U = — 1, 
corresponding to the tip of the lobe. 



energies. It is a quantity of primary importance in cold 
atom experiments, as it is directly accessible by imaging 
the cloud after expansion^i^iH. We have considered a 
two-dimensional system, which allows for the existence 
of Bose-Einstein condensation with long range order in 
the SF regime. 



The 2D momentum distribution is shown in Fig. [8] for 
different values of the parameter t/U. We have checked 
that we can reproduce the momentum distribution in the 
two extreme cases of deep MI (see Eq. ©) and weakly 
interacting SF (see Eq. (fT^ ). 

In the Mott phase with non vanishing tunneling, we 
find that the momentum distribution presents a modu- 
lation showing up as interference peaks in the expansion 
pictures^^i^. When the QPT is reached, a strong peak 
develops at q — corresponding to the condensate. How- 
ever, in a SF close to the phase transition the background 
momentum distribution (at g 0) is large indicating a 
strong depletion from the condensate due to interactions. 

The momentum distribution obtained with the RPA 
method happens to be not correctly normalized to the 
total number of atoms. We attribute this feature to the 
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FIG. 8: (Color online) Momentum distribution n{q) for a 2D 
system for /x/C/ = 0.5 and (a) t/U = 0.05 deep in the MI; (b) 
t/U = 0.15 in the MI but closer to the QPT; (c) t/U = 0.175 
in the SF phase close to the QPT; (d) t/U = 0.25 in the SF 
phase but further from the QPT. 



fact that fluctuations are not self-consistently included in 
the ground state (see discussion in App. [C|). 

Starting from the momentum distribution, one can ob- 
tain the single particle density matrix p(r, r'), which con- 
tains direct information about the spatial correlations 
present in the system. In the SF phase, the system 
is characterized by off-diagonal long range order, and 
at long distances, the single particle density matrix ap- 
proaches a constant value equal to the square of the order 
parameter, or the condensate density no 



p(r,r') = (a^(r)a(0))^¥''-^o, 



(17) 



which is non-zero in a SF. This quantity is linked through 
Fourier transform to the i5-function at q = which ap- 
pears in the momentum distribution. 



The single particle density matrix in Fig. [9] shows a 
marked transition from the MI phase to the SF, and cor- 
responds to a change from an exponential decay of the 
correlations to a (quasi) long range order. In the MI, the 
single particle density matrix shows that the correlations 
decay over a finite length scale, which decreases by de- 
creasing the tunneling and moving deeper into the Mott 
lobe. From those results, we have direct access to the 
length scale of the correlations in the insulating regime. 

The condensate fraction can be in principle extracted 
from the momentum distribution, by subtracting the de- 
pletion X]q/o ""(l) fro™ the total density n, or equiv- 
alently by looking at the asymptotic value at large dis- 
tances of the single-particle density matrix p(r, r'). How- 
ever, unfortunately we find that the present application 
of RPA gives a violation of the total density sum rule 



(b) 



20 

r — Tl 



-20 20 

r — rl 





FIG. 9: (Color online) Density matrix p(r, r') as a function 
of the relative distance r — r' for a 2D system (black line 
for a cut in the center of the trap and blue line along the 
diagonal), (a) t/U = 0.05 deep in the MI showing only nearest 
neighbor correlations; (b) t/U = 0.15 in the MI but closer to 
the QPT showing an increase in the scale of the short range 
correlations; (c) t/U = 0.175 in the SF phase close to the 
QPT showing long range order; (d) t/U — 0.25 in the SF 
phase but further from the QPT. Note that the asymptotic 
value for large r — r' has been subtracted. 



and J2q^o''^i'i) > ^Iso discussion in App. [C]), so 

that neither the momentum distribution nor the single- 
particle density matrix turn out to be useful quantities to 
extract the condensate density. This problem arises be- 
cause within the present theoretical description, we have 
not included the feedback of the collective modes and 
other excitations into the mean field ground state. Pos- 
sible solutions to this problems will be topic of further 
research'^?. 

However, our analysis of the excitation modes and their 
strength allows us to extract the condensate density uq 
directly, using our knowledge of the sound velocity (from 
the slope of the phononic mode at small momenta) , com- 
bined with the knowledge of the strengths of the spectral 
function and the compressibility. Starting from the rela- 
tion 



Ps 



(18) 



and assuming that the condensate and SF densities are 
equal (at least within mean field theory), we get 



no 



'm*dii/dp d/i/dp 



(19) 



where Sph is the strength of the phononic modes. We 
have used t = h^/{m*d'^) (with h = l/2n and lattice 
spacing d = 1) to relate the effective mass and the tun- 
neling parameter in the BHM. 
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The limiting behavior of the RHS of Eq. (|T9|) for very 
small q approaches the condensate density which coin- 
cides with the predictions of mean- field theory (|<y5p), as 
shown in Fig. [TUl 




0.4 0.6 
q/7l 

FIG. 10: (Color online) Condensate fraction obtained from 
the strength of the poles of the spectral function, as in 
Eq. (|19p . for the phonon mode at positive (upper green line) 
and negative (lower red line) frequency. As g ^ both func- 
tions approach the condensate fraction no obtained with MFT 
(dashed blue line). In this figure iJ./U = 0.5 and t/U = 0.5. 
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FIG. 11: (Color online) Inhomogeneous system: (a) order 
parameter in the trap: It is zero in the central MI core and 
finite and large in the surrounding SF ring. The other panels 
show G{v, r') for a given value of the energy as a function of 
r' for a fixed r (black dot): (b) r = (—6, 0) in the SF ring; (c) 
r = (—3,0) at the interface between the MI core and the SF 
ring; (d) r = (0, 0) in the center of the MI core. 



VI. INHOMOGENEOUS SYSTEM: OPTICAL 
LATTICES IN AN EXTERNAL TRAPPING 
POTENTIAL 

We extend the RPA formalism to real space and in- 
clude a spatially inhomogeneous potential, which is taken 
into account in the BHM through a site-dependent chem- 
ical potential. The self consistent MF solution produces 
alternating shells of insulating and superfluid phases 
moving out from the center to the edge of the trap^ii^i^. 

In the inhomogeneous system, the derivation of the 
equation for the Green's function is the same outlined 
in App. [Al with the essential caution that the on-site 
energies and the tunneling coefficients T^^^^, depends 
on position and must be calculated respectively for each 
site and for each pair of neighboring sites. 

In the presence of an external trapping potential, the 
density and order parameter become non-uniform as 
shown in Fig. [TIT a). With our specific choice of param- 
eters, one finds a central Mott core at density n = 1, 
surrounded by a ring of superfluid. The sequence of pan- 
els (b) to (d) shows G(r, r', w) as a function of r' for fixed 
r and oj, which roughly speaking represent the effect of 
perturbing the system at different points: perturbations 
in the SF regions produce a large effect all along the SF 
ring. As the perturbation moves to regions with lower or- 
der parameter near the SF-Mott interface its effect gets 
reduced and finally perturbations in Mott-like regions de- 
crease exponentially and produce negligible effects. 



The results shown in Fig. [TT] are obtained for a given 
value of the energy w. At different energies the structure 
remains the same, but the period of the oscillations along 
the ring changes. By integrating G(r, r',w) over uj, one 
gets the equal-time correlation function G(r, r'). This 
quantity would maintain the ring-like behavior shown 
Fig. [Tl] and hence be qualitatively similar to the one cal- 
culated by Wessel and collaborators'*^. 

Before searching for quantitative results in the non uni- 
form system, one should ponder on the consequences of 
the problem in the normalization of the momentum dis- 
tribution discussed in Sect. |Vl which might affect the 
extension of the RPA method to trapped systems. While 
this extension is technically simple, although it might be- 
come computationally quite expensive, its validity is to 
be questioned due to the co-existence in the same sys- 
tem of different phases (MI and SF), where, as we have 
explained above, RPA introduces different normalization 
factors. For this reason, while we believe one can get 
some insights on the correlations in the system, those 
pieces of information can be trusted only at the qualita- 
tive level. 



VII. CONCLUDING REMARKS: RPA VS 
MEAN-FIELD 

In this work, we have studied the excitations and the 
spatial correlations of the BHM in the RPA approxima- 
tion. It is interesting to compare the information that is 
obtained between the mean field approximation and the 
RPA which includes a certain class of fluctuations: 

(1) The prediction of the SF to MI transition in the 
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[/i/?7, t/U] plane calculated from the vanishing of the or- 
der parameter tp at the mean-field level and by the dis- 
appearance of the gapless excitation mode within RPA 
lead to exactly the same result for the boundaries of the 
insulating lobes. 

(2) At the mean field level, one can extract informa- 
tion about the sound velocity using the compressibility- 
effective mass relation, while at the RPA level the sound 
velocity is given by the slope of the phonon mode. The 
two methods again lead to exactly the same result. 

(3) The condensate fraction is no = |<^P at the mean- 
field level, whereas in the RPA treatment, it is extracted 
from the small q behavior of the spectral function and us- 
ing the mean-field compressibility. Again the two meth- 
ods lead to exactly the same result. 

If the fluctuations around the mean field state were in- 
cluded self-consistently, they would renormalize the en- 
ergy and order parameter. We then expect the SF state 
to be susceptible to fluctuations and the critical t/U to 
be shifted to a higher value than the one obtained by the 
MF theory 

The RPA method further gives information which is 
not included in the mean-field treatment. These include 
(i) the excitation spectra; (ii) their strengths; (iii) the 
existence of new gapped modes in the strongly interacting 
SF phase; (iv) the momentum distribution and spatial 
correlations in the system. 

As an open question, we are left with the role of quan- 
tum fluctuations in the vicinity of the QPT, which as we 
discussed may have additional effects on physical quan- 
tities like the e.g. momentum distribution or the con- 
densate fraction. As explained in App. [Cl it is to be 
expected that in our approach the momentum distribu- 
tion is not normalized. In the Mott limit, a clear devia- 
tion from the normalization to the total number of atoms 
can be calculated analytically. Analogously, in the Bo- 
goliubov regime, exactly recovered in RPA in the dilute 
limit, where one assumes n = the normalization is 
larger that the total number of atoms once the depletion 
(given by the integral over all momenta different from 
zero of is added. However, while in the deep MI 

and SF regimes, the change in the normalization is just 
a small perturbation, close to the phase transition it is 
a striking effect. We attribute this to the fact that we 
perform RPA on the mean-field ground state, without 
taking the effect of RPA self-consistently into account. 
To this same reason, we attribute the existence of pre- 
dictions (1) to (3) which are equal in the mean-field and 
RPA treatments. 



APPENDIX A: GREEN'S FUNCTION 
FORMALISMS IN THE RPA APPROXIMATION 

In this appendix we recall the main steps of the deriva- 
tion of the Green's function formalism in the random 
phase approximation (RPA)^i21. 

RPA includes some fluctuations around the mean-field 
solution, which allows us to describe the excitations of 
the system. However, as explained in Sect. IVIIl these 
fiuctuations are not included self consistently, allowing 
feedback into the mean field ground state. Also ignored 
are quantum fiuctuations of the order parameter which 
are especially important close to the QPT (see also dis- 
cussion in App. [Cjl . 

Mean-field decoupling: 



Substituting a = ip + a, = Lp 
obtain without any approximation 



a\ in Eq. (P), we 



^Uiijii - 1) - ^iUi - tip{a\ + ai) + tip^ 



t 



(Al) 



The Hamiltonian H is thus rewritten as a sum of on-site 
Hamiltonians Hf"^^ (indicated by the term in the square 
bracket which includes hopping at the mean field level), 
plus an inter-site hopping term, which is assumed to be 
small. 

Random Phase Approximation: 

In the basis given by the complete and orthonormal 
set of on-site eigenstates \ia) of the on-site Hamiltonians 
H^^ , the Hamiltonian in Eq. (|Aip takes the form 



— E^a^"a~2z E '^aa't3l3'^aa'L0p,, (A2) 



{ij)aa'Pf3' 



where we have defined 



\ia){ia'\, 

{ia\a\\ia'){jP\a,\jP') 



{ia\a,\ia'){]l3\ii]\jl3' 



(A3) 



For any pair of single particle operators A and B, the re- 
tarded (77 = -f 1) or advanced (77 = —1) Green's function, 
defined as 
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a,a(r) = -iTiO{r,T){A{T)B{Q) - S(0)A(r)), (A4) 
can be written in the on-site eigenbasis as 



G(r)= {^a\A\^a'){JmJP)G'Lf^p'ir), (A5) 
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where G'^^,p^,{T) = ((LL,(t); L^^^,)). 

In the energy domain, the Green's function is defined 



G{u;) ^ {{A-B))^^ = I ^G(r),,,e-±^ (A6) 



where 



A relation analogous to Eq. (|A5|1 also holds for the energy 
resolved Green's functions. 

In the uniform system at vanishing tunneling (deep 
Mott regime), H = ^ e^^La is exactly diagonahzed by 
the on-site eigenstates (Fock basis). The Green's function 
G{uj) can be calculated exactly to be 



n(q,^) 



An = 



A22 
A12 

A21 



An 



E 



Ai2A2ie{q) 
+ 1 - e(q)A22 ' 



LO — A„ 



oj — Aq, 



^12; 



(All) 

(A12) 

(A13) 
(A14) 

(A15) 



1 



1 



(S« 



{En-l — E„ 



(A7) 



with En = -nn + {U/2)n{n - 1). 

When the nearest neighbor tunneling plays a role, the 
commutator with the transition operators L^^^, with the 
tunneling part of the Hamiltonian produces a coupling 
to three operators Green's functions. Following the pre- 
scriptions of RPA of replacing the average of a product 
with product of averages, one obtains 



{E-el^,+el 



(A8) 



ZTT 



{fc)i77' 



At zero temperature results, are equal to 1 for the 

ground state and vanish otherwise. 

For nearest-neighbor hopping and for a uniform sys- 
tem, where e^, and f^'f^^y are site- independent 
(see Eq. (jABI) ). the same equation in momentum space 
takes the form 



Eq) Gaa'P0'{E, q) = 
-{{Laa) — {La'a'))Sa0'Sa'(3 + 



{E-e^., +ea)Gao.>p0'{E,q)= (A9) 

1 

+ {{Laa) — (-^g'g')) Tg, a^^' G ^^i p fji {E , q). 

77' 

where e(q) — — {2t/z) cos{qi), with i running over the 
dimensionality of the system and z being the number 
of nearest neighbors. In practice, for each value of the 
energy E and momentum q, the solution of Eq. (jA9p 
amounts to inverting a 2(iVs — 1) x 2(iVs — 1) matrix, where 
Ng is the dimension of the number state basis considered. 
The solution can also be found analytically to bei^ 



G(q,c^) 



1 



n(q,c^) 



27r 1 - e(q)n(q,w)^ 



(AlO) 



with |0) indicating the ground state, j/^q — {a\a^\0) (and 
analogously the other terms), and — Ea — Eq. This 
equation can be easily evaluated numerically once the 
mean-field ground state wave function is known. This 
can be done analytically in the MI phase (see App. [C|) . 
but has to be done numerically in the SF phase. 

In principle, it is possible to apply this formalisms also 
in the non uniform case. One should then start from 
Eq. (|A8p considering that both and T^^^^, generally 
depend on position. Hence, in that case, for each value 
of the energy the solution amounts to inverting a 2(A'^s — 
1)N X 2{Ns — 1)N matrix, where TV is the number of 
lattice wells considered. 



APPENDIX B: BOGOLIUBOV THEORY FOR 
THE BHM 



In this appendix we will present the details of the 
Bogoliubov treatment for the BHM. The results are ex- 
pected to be valid in the weakly interacting SF regime, 
and were compared to the results of the RPA approxi- 
mation in Sect. IIVI 



We start from the BHM 



I i (ij) 



(Bl) 

and define, as before, the fluctuation operators subtract- 
ing from the operators a and their mean-value 



~ t f * 
a = a — ((9, ~ — ip . 



For a uniform system, one gets 



(B2) 
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i 

i 



(B3) 



3rd + 4th order in 5 and 



To minimize the energy one has to set to zero the first 
order, leading to the discretized version of the Gross- 
Pitaevskii equation for the uniform system 



U\ip\^ - fi - t = 



^GP 



u 



(B4) 



The quantity \ip\'^ is the density in the uniform system, 
and is related to the chemical potential. Conversely, for 
a given density the chemical potential is given by 



^l^u\ipf-t, 



(B5) 



namely the kinetic energy plus interaction energy. 
Strictly speaking the kinetic energy of the condensate is 
zero, and t is just a shift due to the choice of the zero of 
energy. This approach does not include any phase tran- 
sition, because the SF fraction is always equal to the 
total density. 

To diagonalize the second order terms in the Hamilto- 
nian {H2), we perform a transformation to momentum 
space 



q 

-E 



(B6) 
(B7) 



so that, finally it reads 



^2 =-^^[2C/|^p-M + e(q)] + 
q 

+ i^[(2C/|(^p-/. + e(q))(a^5, 



(B8) 



*2~ 



where, as before, e(q) = — (2t/z) cos(gi). Then, we 
apply the Bogoliubov transformation which diagonalizes 

^2 



Uq6q 



-q^'-q 



■ q^-q + Wq^q 



(B9) 



with the additional condition — |f-qP = 1 to pre- 
serve the commutation relations. This is equivalent to 
the Bogoliubov equations 



(£q - ftWq)Uq + [/(/J^Wq = 0, 



(£q -I- fiiJq)Uq + U ^p* U. 



0, 



where £q = U\ip\^ + {Mjz) sm\q,d/2). 

The solution for the Bogoliubov spectrum and the Bo- 
goliubov amplitudes is given by 



Ua+V„ = 




(4Vz)E.sin^fa/2) 

hWq ' 
hujq 

(4i/z)E.sin2(9./2) 



(BIO) 
(Bll) 

.(B12) 



For q ^ 0, the spectrum shows a linear behavior in 
|q|, huq « \q\^y{2t/z)U\^p\'^ with sound velocity c = 

In the Bogoliubov theory the momentum distribution 
is given by 



n(q) = |(y5p5q,0 + b-ql^ 



(B13) 



It is evident that starting from the fact that l^?^ equals 
the total density in the uniform system, the integral of 
the momentum distribution Jn(q)(iq will exceed this 
value by the depletion — J \v{q)\'^dq. In the limit 
of validity of the Bogoliubov approach, this quantity is 
very small. 



APPENDIX C: ANALYTIC CALCULATION OF 
THE MOMENTUM DISTRIBUTION IN THE 
MOTT REGIME IN THE RPA APPROXIMATION 

In the RPA approximation, the momentum distribu- 
tion in the MI regime can be calculated analytically. In 
this appendix, we will present the results for 1,2, and 3D 
systems, with the aim of pointing out that close to the 
phase transition quantum fluctuation play a major role 
and are not correctly taken into account by RPA. 

In the MI regime, the RPA Green's function in 
Eq. (jAlOp can be written as 
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GM/(q,w) 



Aiiiuj) 



27rl-e(q)Aii(w)' 



(CI) 



where An has been defined in Eq. (jAl2|) . A12 = A21 =0 
in the MI, and e(q) = —{2t/z) J^i cos(qi) as before. 

In particular for the MI at density n (i.e. where the 
on-site ground state is |0) ~ n), we get 



An{uj) 



n + 1 



^ — {En+l — En) UJ + {En-1 — E„ 



(C2) 



(with En — ~nfi + {U/2)n{n — 1)) which obviously co- 
incides a part from a factor 1 /27r with the Green's func- 
tion for the MI at zero tunneling previously introduced 

in Eq. ,— ^ 

By using definitions in Eqs. (jCllC2p above and defining 
A± = En±i — En, it is straightforward to see that the 
Green's function takes the form 



Gmi — 
1 



(C3) 



{n + l){uj + A_) - n{Lu ~ A+) 



27r (cj - A+)(cj + A_) - e(q)[nw + {n + 2)A_] ' 

whose poles uj± can be calculated analytically and have a 
momentum dependence due to the kinetic energy e(q) in 
Eq. (jCip . Hence, the Green's function close to the poles 
(lu w <-u±), is approximately 



G 



MI 



1 {n + l){uj± + A_) - n{uj± - A+) 

27r {uj± — LUzf)(LU ~ LU±) 

Consequently, the momentum distribution reads 



(C4) 



n(q) = -(" + ^^("-+^-)-"("--^+). (C5) 



■UJ+ 



In the deep MI {t = 0), uj± = ±A± and the momen- 
tum distribution is simply given by n(q) — n, which in- 
tegrated over the allowed momenta (equal to the total 



number of wells) gives the total number of atoms. At fi- 
nite tunneling, the kinetic energy e(q) gives a modulation 
to the momentum distribution. 

To find the normalization of n(q), we integrate numer- 
ically the analytic expression in Eq. (jCSp for different di- 
mensions. The result is shown in Fig. [T^] for the MI with 
n — 1, and clearly shows an increase of the normaliza- 
tion toward the phase transition. We attribute this effect 
to the fact that quantum fluctuations are not completely 
taken into account in this approach. This interpretation 
is confirmed by the fact that this effect diminishes by 
increasing the dimensionality of the system. 
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FIG. 12: Normalization of n{q) for jj,/U — 0.5 as a function 
of t/(7 in dimensions ID (blue dashed-dotted line), 2D (green 
dashed line) and 3D (red full line). The thin vertical line 
indicates the phase transition. 



It is interesting to note that while the normalization of 
the momentum distribution is violated, the sum-rule is 
perfectly satisfied. The strengths at positive and negative 
energy are respectively given by 



such that iS"- 



(n + l)(w± + A_)-n(wj 



UJ± — UJ^ 



(C6) 



S- = 1. 
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